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ABSTRACT 

Modeling the late inspiral and merger of supermassive black holes is central to understanding accretion 
processes and the conditions under which electromagnetic emission accompanies gravitational waves. We use 
fully general relativistic, hydrodynamics simulations to investigate how electromagnetic signatures correlate 
with black hole spins, mass ratios, and the gaseous environment in this final phase of binary evolution. In all 
scenarios, we find some form of characteristic electromagnetic variability whose detailed pattern depends on 
the spins and binary mass ratios. Binaries in hot accretion flows exhibit a flare followed by a sudden drop in 
luminosity associated with the plunge and merger, as well as quasi-periodic oscillations correlated with the 
gravitational waves during the inspiral. Conversely, circumbinary disk systems are characterized by a low 
luminosity of variable emission, suggesting challenging prospects for their detection. 

Subject headings: accretion, accretion disks — black hole physics — gravitational waves 



1. INTRODUCTION 

Gravitational waves (GWs) from the inspiral and coa- 
lescence of supermassive binary black holes (BBHs) con- 
tain detailed information about the black hole (BH) masses, 
their spins, and the orbital characteristics of the binary. 
Space-based GW observations will provide measurements of 
these quantities with a lev el of accuracy rarely attained by 
astronomical observations dTrias & Sintes|[2008t iKlein et alj 
|2009). Electromagnetic (EM) observations of these cata- 
clysmic events, on the other hand, can provide an observa- 
tional link between the merging BHs and their host galaxies 
by shedding light on the environment surrounding the holes 
and, in particular, by giving us insight about accretion and 
feedback processes. Our work investigates the circumstances 
under which detectable EM emission accompanies the GW 
signal from supermassive BBH mergers in astrophysical en- 
vironments. 

Our work complements non-relativistic hydrodynamic sim- 
ulations that follow the inspiral of a supermassive BH pair in 
the aftermath of galactic merger down to a scale of O.Olpc. 
(see|Q)lpi et al. 2007, for a review). These studies have pro- 
vided clues about BH interactions with the surrounding stars 
and gas that lead to the formation of gravitationally bound bi- 
naries. They suggest that characteristic EM signatures may 
be associated with supermassive BBHs in this stage of evo- 
lution ( Armitage & Nataraian 2002; Milosavlievic & Phinney 
120051: iDotti et alj 120061: iBogdanovic et alj |2008, for exam- 
ple). Modeling at even smaller scales, where the orbital dy- 
namics of the binary is determined by gravitational radiation, 
requires a post-Newtonian treatment of gravity and, in the 
last few tens of orbits and merger, the framework of gen- 
eral relativity. Examples of such works include recent nu- 
merical relativ ity studies of coalescin g bina ries surrounded by 
test particles (Ivan Meter et all I20T0I) . gas jBode et al.ll2O10i: 
Farris et al. 120 lOI). or electrom agnetic fields jPalenzuela et alj 
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Previously, in iBode et akl (I2010I) . we carried out the first 
fully general relativistic, hydrodynamics study of the late in- 
spiral and merger of binaries with equal-mass and parallel- 
spin BHs in a hot accretion flow. This work together with 
iFarris et aTl (120101) established that plunge and merger of 
"symmetric" binaries in hot accretion flows are characterized 
by an EM flare followed by a sudden drop-off in luminos- 
ity. Another of its findings is that the BBH orbital motion 
imprints a quasi-periodic signal in the light curve that is an 
EM equivalent of a GW chirp. However, it remained to be 
answered whether such signals are present for more generic 
binary configurations and in cooler accretion flows such as 
accretion disks. The present work investigates the effect of 
varying BH spins and mass ratios and considers astrophysi- 
cally relevant environments that include circumbinary disks 
and hot, radiatively inefficient accretion flows. The hot flow 
and circumbinary disk effectively bracket a range of physical 
scenarios for a BBH environment characterized by the bal- 
ance of heating and cooling processes in the gas. If cooling 
is more efficient, the gas settles into a rotationally supported 
accretion disk whereas when heating dominates it gives rise 
to a hot, tenuous, and geometrically thick accretion flow. 

We consider binaries with BH masses m,- {q = 1112/1111 < 1) 
and dimensionless spin parameters a)/m,. The initial binary 
is on a quasi-circular orbit at a separation of 8M (in geo- 
metrical units), with M = m\ +1112 the total mass of the bi- 
nary. We compute our results using M as a natural unit and 
present results scaled to a total mass of lO'' Mq; that is, M = 
L48 X 10'-cmM7 =49.4sM7 withM? =M/1O^M0. Since in 
all scenarios the mass of the gas in the vicinity of the binary 
is many orders of magnitude lower than M, the BBH dynam- 
ics is indistinguishable from the vacuum case. We use the 
computatio nal met hodology and infra st ructure described i n 
Bode et all (12 010); Ansorg et al.' (2004); Baiotti et al. (2005*); 



CactusI (12010); Einstein Toolkit (■2 (310)j_ Hu sa et al. (2006); 
Schnetter etal] (^2004^ ; iTh'ornburgI (|2004 . Our simulations 



do not capture radiative transport, magnetic fields nor physi- 
cal viscosity. Our computational grid had an outer boundary 
at ^ 320M, with 9(10) levels of mesh refinement for the equal 
(unequal) BH mass binaries. Resolution at the finest level was 
M/67(M/76), with subsequent levels increased by a factor of 
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two. The 5 (7) finest levels had 42 (36 ) grid points, while re- 
maining had 84^(72^) grid points. A subset of runs at higher 
resolution of M/ 100 showed that our results did not change 
significantly with resolution to affect our conclusions. 

2. HOT ACCRETION FLOW 

In this scenario, the BBH environment is assumed to 
have the properties of a radiatively inefficient accretion flow 
(RIAF). In RIAFs, most of the energy generated by accre- 
tion and turbulent stresses is stored as thermal energy in the 
gas, and the accretion flow is hot a nd geometrically thick 
(llchimarul [19771: iNaravan &^ 1 19941) . The electron and ion 
plasmas in RIAFs can form a two-temperature flow in which 
the thermal energy is stored in the ion plasma while the elec- 
tron plasma cools more efficiently (i.e., Tp > T^). In such a 
case, the temperature of the plasma is represented by the ion 
temperature, while the characteristics of emitted radiation de- 
pend on the properties of electrons. The temperature ceiling 
reached by the ion plasma will be determined by cooling pro- 
cesses such as the thermal bremsstrahlung, synchrotron, and 
inverse Compton emission, as well as the electron-positron 
pair production and the pion decay resulting from energetic 
proton-proton collisions. Which process dominates the en- 
ergy loss of the plasma sensitively depends on its density, 
temperature, and magnetic field strength, as well as the ef- 
ficiency of coupling between ions and electrons. The latter 
process determines the rate with which energy can be trans- 
ferred from hot ions to electrons, and consequently the ratio of 
their temperatures, e = Te/Tp. While modeling of the radiative 
cooling and treatment of the two-temperature plasma flow is 
beyond the capabilities of our code at this time, we capture the 
effect of different thermal properties of the plasma by inves- 
tigating several initial ion plasma temperatures in our simula- 
tions, Tp = {10"', 10", 10'^}K. In order to evaluate the emis- 
sion properties of these flows, we make a simplistic assump- 
tion that e = 10"^ everywhere in the accretion flow. This is 
an idealization as T^/Tp is expected to vary in both space and 
time and can have a range of values between ^ 10"^ and 0.1 
dep ending on the domin ant plasma processes (see, for exam- 
ple, [Sh^maeral]|2007l)- Our choice of e is however conser- 
vative because it caps and consequently the bremsstrahlung 
luminosity at lower values (see § 12.2b . The properties of the 
binary and initial Tp of the environment are listed in top part 
of Table [U 

2.1. Initial Conditions 

The gas around the binary is initialized with uniform den- 
sity and temperature and modeled with a polytropic index 7 = 
5/3, adequate for the ion plasma which is non-relativistic at 
temperatures Tp < lO'^K {kTp < 100 MeV). The gas initially 
has zero linear and angular momentum. The latter condition 
is consistent with the expectations bas ed on the s elf-sim ilar 
solutions for RIAFs described by ,Naravan & Yil (1 1 9941) for 
7 = 5/3 gas. This initial configuration is first evolved on 
the static BBH spacetime for the duration of 32M and subse- 
quently on a dynamic BBH spacetime for 160MQ In compar- 
ison, the relaxation time of the gas within the Bondi radius of 
the gravitational influence of the binary (freiax = Rb/cs, where 
c, is the speed of sound) is 18M, 530M, and 1.7 x lO'^M for 
10'^ K, 10" K and 10'" K gas, respectively. Longer relax- 

To ensure that this initial phase of relaxation does not introduce spurious 
transients to our analysis, we omit it and only report the subsequent evolution. 
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TABLE 1 

Parameters for the scenarios discussed for hot accretion 
flows (top half) and circumbinary disks (bottom half). the 
central columns describe the the black hole parameters of 
the system while the last column designates the primary 
parameter for the initial gas configuration. the asterisk 
next to dc 1 refers to retrograde circumbinary disk 
rotation. 

ation time scales for lO'^K and 10"K gas imply that at the 
beginning of the simulation the gas is still settling into a quasi- 
hydrostatic equilibrium and flowing toward the center of the 
potential. We verified for all simulations that relaxation is a 
gradual process which can be reproduced in simulations of 
accretion flows with the same properties surrounding a single 
black hole with mass equal to that of the BBH. We use the fact 
that the relaxation process does not introduce rapid transients 
(i.e., on the BBH orbital time scale) in the light curve of either 
system to remove a smooth secular modulation in luminosity 
amplitude caused by relaxation by dividing a light curve cal- 
culated for the BBH system by its single BH equivalent. Us- 
ing this ansatz we disentangle the effects of the gas relaxation 
from variability driven by the orbiting BBH while at the same 
time reducing the computational expense of our simulations 
to about 800M per simulation. This computational gain ulti- 
mately allows us to explore a relatively wide parameter space 
of the BBHs and gas in this work. 

2.2. Properties of the Accretion Flow 

In all models, a pair of denser gas wakes quickly forms be- 
hind the orbiting holes after the beginning of the simulation. 
About two orbits before merger, a high-density, dynamically 
unstable region is also formed inside the binary orbit. Af- 
ter coalescence, this structure is promptly swallowed by the 
newly formed BH. Both the wakes and the high-density cen- 
tral region are associated with strong shocks, driven by the 
orbiting BBH, and their evolution gives rise to the charac- 
teristic variability of the emitted light. This is illustrated in 
the top panels of Figure [T] which show bremsstrahlung lumi- 
nosity as a function of time, integrated within the sphere of 
radius R = 2QM and normalized by the light curve for a single 
BH with equivalent mass. We verified that variability in the 
accretion flow associated with the BBH evolution is enclosed 
within this volume by varying our choice of R and confirmed 
that the central region contributes the dominant component of 
the luminosity to the light curve of the binary. 

The top-left panel in Figure [T] shows a broad peak in lu- 
minosity whose growth coincides with the formation of the 
shocked region within the binary orbit. As the binary shrinks, 
the brightness of this region increases until merger at which 
point the final BH swallows the shocked gas and the lumi- 
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Fig. 1. — Top: bremsstrahlung luminosity in hot accretion flows for different BBH configurations (left) and accretion flow temperatures (right) normalized to 
the luminosity of a single BH with corresponding mass. Bottom: ratio of beamed to unbeamed luminosity L^b for different BBH configurations (left) and different 
accretion flow temperatures (right) for an observer viewing the binary edge-on. Coalescence occurs at f = 0. 



nosity precipitously drops off. The top-left panel of Figure [T] 
compares the light curves calculated for different BBH con- 
figurations in an environment with initial temperature Tp « 
10'^ K. The q = 1/2 case (CH2) exhibits a lower and nar- 
rower luminosity peak relative to the q = I case (CHI). This 
is because orbital torques from an unequal-mass binary are 
less efficient at driving shocks in the gas, resulting in a less 
luminous and delayed emission from this region. In q= 1/2 
systems with generic spin orientations (CH3 and CH4), the 
luminosity peaks at even lower values, at ^ 80% of the am- 
plitude achieved by the q = I binary. This is a consequence 
of the orbital precession present in binaries with misaligned 
spins, which further inhibits the formation of a stable shock 
region between the holes. The gradual rise and sudden drop- 
off in luminosity however seem to be a generic feature of all 
modeled light curves, regardless of the spin configuration and 
the mass ratios explored in our simulations. Moreover, the 
same feature has been observed for a wide range of initial con- 
ditions employed in previous works ca rried out by our gr oup 
jBode et al.l2010l) and by other authors jFarris et alj201(3l) . in- 
dicating that this is a robust signature of binary systems merg- 
ing in hot accretion flows. 

To understand the dependence of the luminosity on the 
properties of the system, we estimate the bremsstrahlung lu- 
minosity emitted from the Bondi radius of gravitational influ- 
ence, /?B 



of the gas or unbind it from the BBH altogether, thus acting 
to erase the variability imprinted by the binary motion and 
suppress the luminosity. 

The top-right panel in Figure [T] shows the light curves cal- 
culated for parallel-spin q=\ binaries in gas with initial tem- 
peratures Tp = {10'°, 10", 10'^} K. Lower temperature flows 
tend to exhibit a larger luminosity peak and a more dramatic 
drop-off than the hotter flows. This is because in colder flows 
the shocks can be very effectively excited by the merging bi- 
nary. We find that regardless of its initial temperature, the gas 
in the vicinity of the holes is persistently shock-heated to a 
temperature ^ 10'~ K as a consequence of the binary orbital 
motion. Hence, the height of the peaks in the top-right panel 
of Figure[T]reflects the ratio of bremsstrahlung emissivities of 
the gas in the BBH system {Tp w 10^^ K) to that in the single 
BH system {Tp « {10'°, 10", 10'^} K), where bremsstrahlung 

1/2 

emissivity cx Tp ' (1 +4.4e_2 Tp^n). It follows from this simple 
estimate that the relative height of the BBH luminosity peak 
is 54, 12, and 1, respectively, consistent within a factor of few 
with the values that we calculate from simulations. 

If the hot accretion flow is threaded by a strong magnetic 
field, a significant fraction of its luminosity could be emitted 
in the form of synchrotron radiation which, assuming field 



strengths = 10"* G/3io''' Tp,n r^'^M/'", could reach 



1/2 ,.-1/2 



;6.5M7;-J2, 



: 4 X lO'^'^ergs 



(2) 



i-brem~6.7 X lO'^^'ergs 
X {\+AAe-2Tp,n] 



-11/2 .^-1/2 
-2 ■'p,12 

5,4 'T '"7 I 



(1) 



where e = 10 e-2, Tp = 10 K Tp u, tj = ctt pRb is the opti- 
cal depth to Thomson scattering within the Bondi sphere, aj 
is the cross section to Thomson scattering, and p is the gas 
density. Subscript "5.4" indicates that the expression in the 
brackets is normalized to this value. Note that Eq. ([TJ implies 
a maximum bremsstrahlung luminosity that can be reached by 
an accretion flow as long as its optical depth tj <l (we con- 
sider Thomson scattering to be the dominant source of opac- 
ity in this regime). Flow with a larger optical depth would be 
subject to radiation pressure which could alter the kinematics 



where /3 = Snpth/B^ = lO/3\o is the ratio of thermal to mag- 
netic pressure in the gas, expecte d to reach values of 1 - 10 
in the central regions of RIAFs dciillolil). The presence 
of the softer photons supplied in situ by synchrotron and 
bremsstrahlung emission would also give rise to the inverse 
Compton radiation of similar magnitude: 



2L,,f,T-i2TTM^, 



(3) 



where relativistic factors have been evaluated for v/c ~ 0.3, 
appropriate for the BBH regime close to the merger, and Ljoft 
is the luminosity of soft radiation. 

Where the high energy tail of protons reaches the thresh- 
old of kTp = 100 MeV an additional high energy pro- 
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cess contributes to the radiative cooling: proton-proton 
collisions result in copious pion production, followed by 
pion decay to two 7-ra y photons, p + p — >■ p + p + 7r° — » 
p + p + 2-r Jmhlbacka e taP [1971 lEilek & KafatosI [1981 
Cploietal. 1986; Mahadevan et al.l Il997t l Oka & Manmoto' 
2003=: .Bhattacharvva et al., ,2006.) . Following IColpi et al.l 
(119861) . who calculated the 7-ray emission from the p-p colli- 
sions of a thermal distribution of protons in vicinity of a single 
Kerr black hole, we estimate 

Lpp « 2- 13 X 10^^ ergs-' Tp,n rf Mj , (4) 

where the two extreme values correspond to a static and max- 
imally rotating black hole, respectively. We expect the lumi- 
nosity in the BBH system to be closer to the higher value be- 
cause the gas in the rotating and dynamic spacetime of the pair 
of orbiting BHs is very efficiently shock-heated to 100 MeV. 
The emission of 7-rays due to the pion decay is strongly sup- 
pressed in the limit T7 > 1 due to the increased cross section of 
7-ray photons to electron-positron pair production, as well as 
the increased coupling between electron and proton plasma, 
which lo wers T„ below the energy threshold for the pion pro- 
duction JColpi et al.lll986l) . In calculating luminosities in this 
section, we assumed the gas to be optically thin within the 
Bondi radius, which sets an upper limit on the gas density of 
the hot accretion flow, and thus Lbrem, ^^syn. Lie, and Lpp. 

The spectral energy distribution of these sources would be 
similar to a group of low l uminosity-AGN to w hich RIAF 
models have been applied jNemmen et al.l [20 06[ for exam- 
ple). Spectral bands where these emission mechanisms are 
expected to peak in the reference frame of the binary are 
submillimeter (synchrotron), UV/X-ray (inverse Compton), 
~ 1 MeV 7-ray (bremsstrahlung and inverse Compton), and 
^ 20 MeV (pion decay). Additional components that we do 
not model in this work but could also arise and in principle 
overtake the emission from the hot gas are the wide-band non- 
thermal synchrotron emission, if active and persistent jets are 
present in the system, as well as the optical/UV emission as- 
sociated with the accretion disk that may encompass the hot 
flow at larger radii |HoH2005t) . 

The bottom panels in Figure [T] focus on oscillations in lu- 
minosity due to relativistic beaming and Doppler boosting of 
light emitted by the gas wakes. We account for these effects 
by multiplying the broadband emissivity of the gas (i.e., the 
luminosity per unit volume) by a factor (W {I - 13 cos{d)))~'^ 
and integrate over the volume to obtain the bolometric lumi- 
nosity. Here, W is the local Lorentz factor and /3cos(6') is 
the projection of the local 3 -velocity to the line of sight of 
the observer. We do not account for the bending of light and 
gravitational redshift of photons in the potential well of the 
BBH. The variations in luminosity shown in Figure[T]are cal- 
culated for a distant observer placed in the plane of the bi- 
nary at infinity, an orientation for which the oscillations in the 
light curve are maximized. To emphasize the oscillations, the 
bottom panels show the ratio of beamed to unbeamed light 
curves. The most notable difference among simulated binary 
configurations is that the q = I case (CHI) yields sinusoidal 
oscillations, while the q= 1/2 cases (CH2, CH3, and CH4) 
give rise to double-peaked oscillations. We also find that in 
configurations with arbitrary spin orientation (CH3 and CH4) 
the sinusoidal peaks exhibit the largest degree of asymme- 
try. This is because the binary with parallel spins (CH2) stirs 
the surrounding gas more uniformly than binaries with mis- 
ahgned spins (CH3 and CH4), which precess about the origi- 



nal orbital plane. 

The bottom-right panel illustrates the dependence of the os- 
cillations on temperature. The most prominent oscillations 
are again associated with the lower temperature flow, which 
is more susceptible to shocking by the binary. While the high 
temperature of the gas in the CHI case prevents formation 
of the strong density gradients, the lower temperature flows 
allow the density wakes to interact and waves of gas to prop- 
agate within the inner region, giving rise to more varied fea- 
tures in the oscillations. In all cases, the oscillations in lu- 
minosity are mirrored by the GWs. Figure |2| shows how the 
double-peaked structure in case CH2 emerges in both the lu- 
minosity and the GWs when the binary is observed edge-on. 
The double peaks in the GWs arise from the superposition of 
I = 2 and / = 3 modes in qj^l binaries. 

Figure [3] shows how oscillations vary as a function of in- 
clination. The light curves have been evaluated for the q = I 
binary for three different temperature cases (CHI, CMl, and 
CLl) and the inclination angle is defined with respect to the 
initial binary orbital plane. Here again we show the ratio of 
beamed to unbeamed light curves. Because the motion of the 
denser gas wakes is tied to the plane of the binary, the oscil- 
lations in all runs disappear with increasing inclination angle. 
The relative drop in the luminosity of beamed light just prior 
to the coalescence is a consequence of de-boosting of light 
caused by the strong radial inflow of the gas toward the black 
holes in this final stage of their evolution. The most notable 
difference among the three runs is that the oscillation curves 
exhibit more structure with decreasing temperature. This can 
be understood because the lower temperature gas is more 
clumpy and more prone to shocks, while in high-temperature 
flows, higher thermal velocity of the gas acts to erase density 
inhomogeneities. Along similar lines, lower temperature gas 
has less pressure support against infall, leading to higher infall 
velocities and thus more significant de-boosting with respect 
to the unbeamed luminosity case. In the precessing binaries, 
the time-varying inclination of the orbital plane imposes an 
additional modulation of the oscillations. Of the situations 
considered in this paper, two generic BBH systems, CH3 and 
CH4, precess due to their misaligned spins. The maximum 
precession angle attained over the entire evolution is at most 
13° with respect to the initial orientation of the orbital plane, 
resulting in only minor modulation in their light curves. 

It is worth noting that accounting for the effects of light 
bending and gravitational redshift would result in the modi- 
fied variability pattern of emitted light, relative to those shown 
in Figure [3] The photons most strongly affected by the gen- 
eral relativistic effects are emitted by the gas wake associ- 
ated with the "background" black hole, as seen from the per- 
spective of an observer This is because these photons travel 
across the deepest part of the potential well of the binary be- 
fore escaping to infinity. However, most of the boosted light 
contributing to the peaks in the oscillations is emitted by the 
foreground wake, associated with the black hole moving to- 
ward the observer and thus with photons which escape from 
the perimeter of the BBH orbit without crossing the deepest 
part of its potential well. These photons will be less affected 
by the general relativistic effects, thus partly justifying our 
simplistic approach in calculation of the light curves. 

3. CIRCUMBINARY DISKS 
3.1. Initial Conditions 
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Fig. 2. — Beamed to unbeamed light curve ratio from case CH2 in Figure[T] 
and corresponding GW in arbitrary units. 
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Fig. 3. — Variation of the bremsstrahlung luminosity oscillations with in- 
clination angle i in accretion flows of initial proton temperature 10'-, lO", 
and lO'" K from top to bottom, respectively. 

In this scenario, we consider circumbinary disks whose in- 
ner edge, the location where the torques from the binary are 
expected to create a low-den sity "hole" in the disk center 
jArtvmowicz & Lubowlll994h . is set at r = 16M. This choice 
is motivated by the requirement that the inner edge of the disk 
"follows" the shrinkage of the binary and remains at twice its 
semimajo r axis until a binary s eparat ion of 8M. In the con- 
text of the lShakura & Sunvaevl(II973h model of a steady-state 
accretion disk, this requirement translates into the value of the 
half-thickness ratio of the disk, h/r, when the rate of viscous 
inflow of the disk, AvIsc = - ^■5{h/r)'^aVK, is set equal to the 
inspiral rate of the binary from gravitational radiation losses 
jPetersll 19641) . Here, Vk is the Keplerian velocity at the inner 
disk edge. This equality is satisfied for li/r = 0.2 (0.4) and 
a viscosity coefficient a = 0.4 (0.1). Note that there are un- 



certainties associated with this estimate since the steady-state 
model does not fully capture the d ynamics of the disk close 
to the BH coalescence. We follow lO'Neill et alj ( 120091) and 
setup gas-pressure-supported disks characterized by a con- 
stant midplane density and entirely azimuthal initial gas ve- 
locities. 

As in the case of hot accretion flows, here we only consider 
BBHs on quasi-circular orbits. This is based on the expec- 
tation that late in the BBH inspiral any initial orbital eccen- 
tricity will be erased due to the emission of gravitational ra- 
diation unless there is a mechanism which can compete with 
it. A mechanism that has been suggested to increase the ec- 
centricity of a BBH orbit is a r esonant interaction of the bi- 
nary with the circumbinary disk ( AFmitage & Nataraianl2005l : 
ICuadra et al]|2009l: iRoedig et alll201 Ih . It is possible to find 
the radius at which the effects of gravitational radiation and 
circumbinary disk are competing by setting e equal for the 
two cases. We a dopt eVisc = 7 ^avisc/2a and the standard ex- 
pression for e'gw jPetersll 19641) . Here, a is the semimajor axis 
of the binary and TZ = a{de/da) is the parameter that relates 
evolution of the eccentricity to that of the semimajor axis. 
lArmitage & NatarajanI (120051) compute TZ from simulations 
and find TZ « -0.1 for hot disks with values of h/r compa- 
rable to those considered in this work. Assuming h/r = 0.2, 
a = 0.4, q=\, and e = 0.8, we find that the two effects balance 
each other at a binary separation of a = 118M. The assumed 
value of the eccentricity follows from a recent finding that 
binary eccentricity driven by the disk t orques alone reache s 
a saturation value in the range 0.6-0.8 (IRoedig et al.ll2011h . 
Following the GW evolution of a BBH on an eccentric orbit 
with e = 0.8 from separation a=118Mto8M, which is a start- 
ing point of our simulations, we find that a maximal residual 
eccentricity is e ?» 0.06. The residual eccentricity only weakly 
depends on the binary mass ratio and has a similar value for 
^ = 1/2 binary. This low residual eccentricity is insufficient 
for the binary to drive shocks in the distant circumbinary disk 
and hence, in the remainder of this paper, we focus our atten- 
tion on the quasi-circular BBHs. 

The parameters of the circumbinary disk models are sum- 
marized in the bottom half of Table[T] All disks are corotating 
with the binary except in model DCl which is retrograde. To 
minimize spurious transients caused by the initial circumbi- 
nary disk relaxing to the dynamic binary potential, the disk is 
"relaxed" for a period of ^ 250M. As before, we do not show 
this initial phase of the simulations as a part of the results re- 
ported here. In the hot accretion flow scenario this relaxation 
effectively eliminates spurious transients. For disks, the re- 
laxation is not as effective and leaves behind residual slow 
oscillations in the bulk of the disk. These do not give rise to 
transients, and have little effect on the properties of the gas in 
the disk hole, where most of the variability is confined. 

3.2. Properties of the Accretion Flow 

In all simulations a plunging binary recedes promptly from 
the inner edge of the disk and as a consequence, the effect 
of the binary on the disk b eyond its inner edge is relatively 
weak. In agreement with lO'Neill et al.l (120091) . we do not 
detect shocks in the body of the disks caused by the binary 
motion nor by the mass loss in GWs associated with the co- 
alescence. In the absence of characteristic variability from 
the body of the disk, we focus our discussion throughout the 
remainder of this section on the EM counterparts emanating 
from the disk hole. 

The two BHs capture a fraction of the gas from the disk 
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Fig. 4. — Density distribution of gas in circumbinary disk models DAI, DA4, and DCl, shown on logarithmic scale, in face-on (top panels) and edge-on 
(bottom panels) view when the binary separation is 2M. White dots mark the positions of the BHs. In panels where only one or no BHs are visible, this happens 
because orbital precession (top middle) or orbital motion (bottom panels) takes BHs outside of the plane of the image. Panels show a region of size 60M. 



inner edge, which flows into the central region to form a low- 
density ambient medium surrounding the binary. Figure |4] 
shows the face-on and edge-on snapshots of the gas density in 
the central region of the DAI, DA4, and DCl cases. Despite 
their differences in mass ratios and spins, the DAI and DA4 
binaries entrain similar amount of gas from the disk. This 
is a combination of the two effects: while ^ = 1 on the one 
hand presents a larger angular momentum barrier for the in- 
flowing gas, it also inspirals at a slower rate than a ^ = 1/2 
binary, which allows it to capture gas over a longer period 
and offset somewhat the primary effect. The q=\ retrograde 
disk (DCl), on the other hand, shows a smaller central disk 
hole, with a higher gas density in it relative to the two pre- 
vious cases. For this case, the binary and the disk interact to 
decrease the angul ar momentum of the gas, as suggested by 
iNixon et"aD (1201 ll) . thus allowing a greater inflow rate. 

To quantify the effect of the binary-disk interaction, we 
calculate the mean density of the gas, ph in the disk hole 
(i.e., R < IQM) as a fraction of the disk midplane density p^, 
namely, /g = ph/pd- For h/r= 0.2 disks, we find /g « 10"^, 
indicating the efficiency of the binary in suppressing the gas 
flow into the central region. For the h/r = 0.4 disk, /g « 10""*, 
owing to the higher thermal speed of the gas which flows in at 
a higher rate. The orientation of spins does not seem to affect 
the gas inflow rate significantly. 

Similar to the hot flow models the luminosity of the emerg- 
ing variable EM signal will be determined by the density of 
the gas in the vicinity of the BBH. In our simulations the 
choice of h/r and a uniquely specify the gas density of the 
disk in hydrostatic equilibrium whose vertical structure is sup- 
ported by the gas thermal pressure. It is worth noting how- 



ever that modeling of the gas-pressure-supported disks in this 
work is a numerical necessity and that the innermost regions 
of realistic accretion disks are ex pected to be predorninantl y 
supported by radiation pressure (IShakura & SunvaevllI973h . 
A key difference between the two classes of disks is that for 
a given a, M, and h/r (or equivalently M) the gas-pressure- 
supported disks are characterized by a higher gas density than 
radiation-pressure-supported disks. This can be intuitively 
understood because with the radiation as a dominant source 
of pressure support, disks need a smaller fraction of the ther- 
mal pressure (and hence, gas) in order to maintain the hy- 
drostatic equilibrium. The implication is that depending on 
the dominant emission mechanism, the luminosity of the gas- 
pressure-supported disks may not trace that of the realistic 
disks for scenarios under consideration. We thus regard the 
density of the disk to be a free parameter, which can be scaled 
to some physical value, and use /g to estimate the fraction of 
the gas captured by the binary. We expect that f„ calculated 
from our simulations is a robust measure, regardless of the as- 
sumed structure of the disk. This is because the accretion of 
the gas from the inner edge of the disk (i.e., its angular mo- 
mentum transport) is likely dominated by the binary torques, 
rather than the radiation pressure or "viscous" effects. This 
point however remains to be confirmed in future simulations 
properly equipped to address it. 

To gain context about the magnitude of variable EM 
signal associated with these systems, let us consider a 
Shakura-Sunyaev model of a radiation-pressure-supported, 
steady-state disk, with h/r = 0.2 in the regime where 
Thomson scattering is th e dominant source of opacity 
( IShakura & Sunvaevlll973h . Using maximum radiative ef- 
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Fig. 5. — Bremsstrahlung luminosity from tlie disk hole (top two panels) in 
units of L|,rem given by Eq. ^5). Total and individual BH accretion rates for 
the case DA2 (bottom panel) in units of Mg given by Eq. ^6). 



ficiency for the Schwarzschild black hole, one finds at the 
disk inner edge a temperature 2.5 x 10^ K and pa ~ 3.5 x 
10~"gcm"''. Since in the central region /g ^ 10"^ and the gas 
is promptly shock-heated to 7), 10" K, the disk hole ac- 
quires the properties of the hot accretion flow discussed in the 
previous section. Assuming that in this case e = Tg/Tp = 0.1, 
the upper limit on the bremsstrahlung luminosity from the 
disk hole and Bondi accretion rate onto the BHs are 



Mb « 2.5 X IO-^Mq yr" Vg,_5 p'^ t;%^ . 
If a strong magnetic field is present, assuming /? = 10, 



(5) 



(6) 



2.5 X lO^^ergs-'/?:,; f's Pa Rio Tp.n M', , (7) 



-1 o-l 



'2 d3 



where e = O.le-i, /g 
R = IQMRw, and T, 



= io"Kr, 



p'a 



= Pd/3.5x 10-"gcm-^ 
p ii. While higher than the 
bremsstrahlung, the synchrotron luminosity is still modest in 
absolute terms forM < 10^ M0 BBH systems. As in the case 
of a hot accretion flow, it can be shown here that the inverse 
Compton luminosity can be of order of Ljyn and Lbrem- In this 
case, however, the shocked gas from the central low-density 
region can more easily (adiabatically) expand into the ambi- 
ent medium. As a consequence, the temperature of the gas 
remains below the 100 MeV threshold for tt" production in 
most of the computational volume and the emission of 7-ray 
photons due to pion decay is expected to be inefficient. 

Figure |5] shows the bremsstrahlung light curves integrated 
within a sphere of radius R = lOM, ensuring contribution only 
from the diffuse gas in the disk hole. The top panel shows 
the q = I cases. The long period variations in the luminos- 
ity in /z/r = 0.2 models (DAI and DCl) correspond precisely 
to the orbital frequency at radius R = lOM and arise due to 
the flow of gas streams on eccentric orbits in and out of the 



integration domain. This is evident in the first and third pan- 
els of Figure in which illustrate that the distribution of gas in 
the vicinity of ^ = 1 binaries is ellipsoidal due to the binary 
torques. The formation of the eccentric streams of matter in 
response to the BBH torques is a salient property of binary- 
circumbinary disk systems and has previousl y been pointed 
out in the context of Newtonian simulations (iHavasaki et al.l 
120071: iMacFadven & Milosavlievicl2008l:ICuadra et al.l2009l) . 
This type of variation is not observed in the "hotter" h/r = 0.4 
disk (DBl). In this case the binary is less efficient at driv- 
ing a resonant motion of the surrounding gas and it finds it- 
self immersed into a relatively uniform density medium. The 
decrease in the overall luminosity in these three cases is con- 
sistent with the scaling Lbrem « implied by Eq. (|5]l. The 
accretion rates reflect a similar behavior, namely, the scaling 
Mb oc /g from Eq. 

The middle panel in Figure |5] shows light curves from the 
^ = 1/2 binaries (DA2 and DA4). The binaries in these cases 
induce a more circularly symmetric distribution of gas within 
the disk hole (see DA4 case in Figure |4]i. The variability or 
oscillations in the light curves are, unlike for q = I binaries, 
associated with the orbiting wakes following the BHs and thus 
follow the BBH orbital motion. Notice also that here again the 
luminosity of the misaligned spins case (DA4) falls below its 
counterpart with parallel spins (DA2). However, unlike the 
case of the hot accretion flows there is no pre-merger flare 
and subsequent drop-off in luminosity at merger. The lack of 
these two signatures in the disk scenario can be attributed to 
the contrasting density distributions in the two scenarios: in 
the hot flow the gas density increases toward the center of the 
binary potential, while in the circumbinary disk case the den- 
sity is lowest at the center of the system, as illustrated in the 
bottom panels of Figure|4] This occurs because the gas is only 
slowly "leaking" from the inner edge of the disk into the cen- 
tral region. The luminosity of a dilute gas powered by shocks 
from the BBH motion just before the merger is overtaken by 
the luminosity of the remainder of the gas in the disk hole. 
Consequently, the luminosity peak associated with the BBH 
merger is invisible in the light curves shown in the top panels 
of Figure |5] 

The bottom panel in Figure|5]shows the accretion rates onto 
the primary and secondary BHs, as well as the total, in model 
DA2. Note that initially the accretion rate onto the lower mass 
secondary BH exceeds that of the primary, because it orbits 
closer to t he inner edge of the disk where it captures and ac - 
cretes gas dArtvmowicz & Lubow|[T996t iGould & Rixll2000l) . 
Later in the inspiral, the primary increasingly captures gas 
bound to the secondary and augments its accretion rate until 
time w -250M before merger, when the accretion rates of the 
primary and secondary are reversed and the primary becomes 
the dominant accretor. 

Our simulations indicate that low-luminosity variability de- 
scribed in this section would be challenging to detect in obser- 
vational campaigns searching for BBHs, except perhaps for 
massive systems with M 1O^M0. The situation is more 
severe for "thinner" disks in which the central regions are ex- 
pected to be even less luminous. In these disks, the radial in- 
flow velocity at the inner edge is oc (/z/ r)^, resulting in the de- 
coupling of the binary from the disk even earlier in the inspi- 
ral. As a consequence, the gas in the disk hole will have lower 
density and thus be even dimmer. The circumbinary disk it- 
self will likely overwhelm the overall emission and shield any 
sign of variability from the disk hole although, strictly speak- 
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ing, the light emitted by the disk is expected to peak in the 
optical and UV-band, while the bremsstrahlung and inverse 
Compton emission from the Tp ~ 10"K gas in the hole is 
expected to peak at ~ lOOkeV. Because of its radiatively in- 
efficient nature the gas in the disk hole would exhibit similar 
emission properties to the hot flow described in § 12.21 



4. CONCLUSIONS 

We studied the EM signatures associated with the late inspi- 
ral and merger of supermassive BBHs with unequal masses, 
different BH spin orientations, and two environments (hot ac- 
cretion flows and circumbinary disks). 

In a hot accretion flow, the plunge and merger of the binary 
give rise to a characteristic flare followed by a su dden drop- 
off. This and earlier hy drodynamic simulations (iBode et al.l 
120101: iFarris et al]l2010l) indicate that the flare is a robust sig- 
nature that arises in all modeled BBH configurations and hot 
accretion flows regardless of the specific setup of initial con- 
ditions. The amplitude of the flare and the drop-off are more 
pronounced in lower temperature flows, where shocking is 
more predominant. We find that the shape and amplitude of 
the flare depend on a mass ratio and spin configuration to 
a lesser degree. A hot flow is also characterized by quasi- 
periodic variability from the beaming of light emitted by the 
gas wakes behind the BHs. In all cases, this EM equivalent 
of a "chirp" is directly mirrored by the GWs (existence of 
such signature has previously been predicted by iKocsis et alj 
120081) . If observed, either signature (flare or variability) would 
be a smoking gun for an impending BBH coalescence in a hot 
accretion flow environment. 

In the circumbinary disk case we modeled geometrically 
thick disks characterized by high radial inflow velocities, 
which can follow the binary until late in the GW-driven inspi- 
ral. We find that a combination of the properties of such disks 
and the prompt decoupling of the binary from their inner edge 
result in an absence of shocks and luminous EM flares from 
the body of the disks. The only region where variable EM 
signatures emerge is the central low-density hole surrounding 
the binary. We find that the gas in this central region also has 
the properties of a radiatively inefficient accretion flow. In 



unequal-mass BBH systems, its luminosity exhibits variabil- 
ity related to the BBH orbital dynamics. In equal-mass sys- 
tems however, the variability traces the dynamics of the gas 
streams plunging from the inner edge of the disk and is not 
correlated with the frequency of GWs. In all cases, the emis- 
sion from the central region is low and based on our models, 
unlikely to be observed for BBH systems < lO^M©. 

In summary, while BBH mergers in both the hot accretion 
flows and circumbinary disks exhibit characteristic EM signa- 
tures, we find that the former is more likely to be sufficiently 
luminous to be observed. To be observable, any variability 
must be more pronounced than the natural variability of "nor- 
mal" AGNs. While the flare signaling the merger seems to 
satisfy this condition, the quasi-periodic oscillations will be 
more challenging to observe. A precursor GW detection by a 
space interferometer could in principle alert EM observatories 
in advance, and hopefully alleviate this challenge. 

An interesting question is whether, in the absence of a GW 
precursor, the flares and variability can be detected in stand- 
alone EM observations. Given that more massive BBH sys- 
tems are also expected to be more luminous, a future serendip- 
itous discovery of > 10*^ BBH coalescence by a burst alert 
mission cannot be excluded. A more systematic search will 
require deep monitoring of the transient sky with multiwave- 
length synoptic sky surveys. While this biases EM searches 
toward the high BH masses, they are complimentary to fu- 
ture GW observations that will likely be optimized to search 
for the lower mass end supermassive BBHs. Both will be re- 
quired in order to eventually understand the properties of the 
BBH population and their role in evolution of galaxies and 
structure in the universe. 
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